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Abstract 

Quantum dynamics simulations can be improved using novel quasiprobability dis- 
tributions based on non-orthogonal hermitian kernel operators. This introduces ar- 
bitrary functions (gauges) into the stochastic equations, which can be used to tailor 
them for improved calculations. A possible application to full quantum dynamic 
simulations of BEC's is presented. 
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One of the oldest problems in quantum physics 
is also conceptually the simplest. How does one 
calculate the quantum dynamical time evolu- 
tion of many-body or strongly interacting sys- 
tems? In this paper, we will treat some recent 
progress towards solving this problem. This uses 
a novel technique of stochastic gauge fields. We 
will focus here on a relatively simple example, 
which allows us to compare numerical results 
with an exact solution. The present results show 
dramatic improvements in sampling error com- 
pared to the previous positive-P [1] distribution 
methods. 

The chief difficulty in many-body quantum dy- 
namics, is that the relevant Hilbert space of 
all but the most trivial cases is typically enor- 
mous. For example, the formation of a small 
Bose-Einstein condensate [2] may easily involve 
N = 1000 atoms with M = 1000 modes, giv- 
ing 10 600 participating quantum states. Simi- 
lar problems also occur in the static calculation 
of many-body ground states and thermal equi- 



librium ensembles. These problems have been 
solved by the use of methods called quantum 
Monte-Carlo techniques [3,4]. It is noteworthy 
that the difficulty of a large Hilbert space is ex- 
actly the same in both the dynamic and static 
calculations. Thus, we have to conclude that di- 
mensionality is not an insuperable barrier. 



The main technique treated here is a class of 
stochastic methods which sample the Hilbert 
space, rather than storing every element of a 
quantum dynamical problem. Provided sam- 
pling errors can be controlled, there is no reason 
why stochastic methods shouldn't be used for 
quantum dynamics, just as they are in QMC 
[3-5] methods used for calculating ground-state 
or thermal equilibrium properties. We present 
methods that are a great improvement on the 
previously used positive-P simulation method 
[1], which is most useful for open systems cou- 
pled to damping reservoirs. In comparison to an 
earlier approach of modifying the noise terms 
dynamically [6] , we focus on methods that allow 
the drift terms responsible for the deterministic 
evolution, to be changed. 



The approach used here is to expand the quan- 
tum density matrix using non-orthogonal co- 
herent state projection operators, together with 
a phase term. This allows a choice of time- 
evolution equations to be made in a way that 
minimizes the phase-oscillations that would 
otherwise occur in a direct path-integral ap- 
proach, while still preventing the phase-space 
oscillations that can occur in the positive-P 
method. 



2 The anharmonic oscillator: a tractable 
model system 



the evolution of the Y-quadrature observable: 
Y(t) = ([a — a)]/(2i)), given an initial coherent 
state. 



3 Hermitian P-distribution 

The positive-P expansion of the density matrix 
p uses a kernel of non-Hermitian coherent-state 
projection operators. Instead, consider a P-like 
distribution with a Hermitian kernel: 

p = J P{a, p, 6, t) A e~ 9 d 2N a d 2N (3 dB, (2) 



A very successful method for time-domain simu- 
lations of damped quantum systems is the posi- 
tive P-representation used in quantum optics. In 
this method, the quantum state is expanded us- 
ing non-orthogonal coherent states. This allows 
multi-boson and multi-mode interacting quan- 
tum systems to be simulated as stochastic pro- 
cesses in the time-domain. These methods have 
been applied to quantum solitons [7], BEC phase 
fluctuations [8] , and to the theory of evaporative 
cooling [9] — where the theory correctly repro- 
duces the formation of a BEC, as observed in 
experiment [2]. 

However, the positive P-representation usually 
has large sampling errors for times after the 
BEC has condensed. This is typical for this 
method, which is most effective for open sys- 
tems coupled to reservoirs. For this reason, the 
remainder of the research presented here has 
been into methods of minimizing the sampling 
error for a very simplified, one-mode version of 
the BEC Hamiltonian: 



with kernel: 



2 V ' 



(1) 



Here a) is the creation operator for a single 
mode of the boson field, with a positive scat- 
tering length constant. The exact solutions 
for some observables can be found directly for 
this simple case, which is clearly very helpful 
while investigating errors. We will focus on 



A = e iy ||a)( / 3|| +h.c, 

e 9 = Tr[A] = 2e Hr cos(# + m). 



(3) 



Here, \ \a) = exp(yjaja|)|0) is an un-normalized 
coherent state, 9 is a real variable representing a 
quantum phase, and n = n r + irii — a - (3*. Any 
state can be represented with a positive hermi- 
tian P-distribution, and expectation values of 
an observable like Y can be calculated accord- 
ing to averages over P. For example, in the one- 
mode case, if the initial condition is a coherent 
state with p — ||ao)(o7 ||, then we expect that: 



?) = (Tr[?A]/Tr[A] 



= Im ( 



traj. 



chq exp 



(e~* - 1) - it/2]) .(4) 



\a \ 2 (e- i! 



4 Stochastic gauges 



Let us now apply the hermitian P-distribution 
to the case of the anharmonic oscillator. The 
master equation is 



(5) 



The next step is to note that there are a num- 
ber of operator identities between terms in the 



Hamiltonian and differential operations on the 
kernel. The ones of interest for this system are 



It is convenient to define 9 = 9 + rii, and to 
introduce the functions G(F) and G(F) : 



(6) G = F + \[n t 



and its adjoint. By using this identity, it is pos- 
sible to transform the operator equation into 
a corresponding Fokker-Planck equation for P. 
First, we change to the more convenient vari- 
ables and ip defined by: 



a = exp 
(3 = exp 



1 -i 
2 

1 - i 



To take advantage of the new phase variable 9, 
consider that the hermitian P-distribution ker- 
nel A also obeys a number of additional differ- 
ential identities in 9. In particular: 



(10) 



G = 



F + \ [rii + n r 



When these are zero the equations are iden- 
tical to those obtained using the positive P- 
distribution. Converting the differential equa- 
tion in P to Ito stochastic equations, we obtain: 



(7) dcp = [n(l - i) - 2G(T + i)\ dt + V2dW, 



dip ■ 



n 



d9 = -2T 



- i) - 2G{T - i)] dt + V2dW, (11) 
G 2 + G 2 ] dt + V2 (GdW - Gdw) . 



The noises dW and dW are random, Gaussian, 
mutually uncorrelated, and uncorrelated for dif- 
ferent times, with variance (dW(t)dW(t)) = dt. 
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A = 0, 5 Anharmonic oscillator with stochastic 
gauges 



E 2 [ ^ + 1 1 A : 



0. 



(8) 



Since these are equal to zero, any multiple of 
them can be added to the master equation with 
no effect, so we have multiplied them by the 
completely arbitrary functions F, F, E 2 , which 
can be dependent on <p,<p* ,ip,ip* ,9 and t. In 
the Fokker-Planck equation formalism, these 
become correspondences for zero. For example, 
defining T = tan(0 + rii) we have: 
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(E 2 P). 



(9) 



Since G and G are completely arbitrary, they 
can be used to tailor the equations to our liking, 
without changing the final results. This is akin 
to what is done with electromagnetic gauges, 
which is why we refer to the G's as itstochastic 
gauges. A suitable gauge, with a free parameter 
H is as follows: 



G 
G 



[rii - n r + \a\ 2 }, 
[rii + rir- |/3| 2 ]. 



(12) 



These correspondences can be added in any 
amount without disturbing the dynamics, as 
long as the boundary terms from partial inte- 
gration vanish [10]. We now wish to convert the 
Fokker-Planck equation to stochastic Langevin 
equations. To do this, the diffusion must be 
positive, and hence we choose: E 2 = F 2 + F 2 . 



The results of simulating the one-mode anhar- 
monic oscillator with this gauge (with two dif- 
ferent values of fi) are shown in Fig. 1, together 
with the positive P results. It can be seen that 
the sampling error in the quadratures has been 
contained and reduced by over twenty orders of 
magnitude! 
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Fig. 1. Expectation value and variance of the 
Y quadrature for the anharmonic oscillator with 
p = |3)(3| at * = 0. The positive P (/x = 0) is 
shown by the dotted line, the hermitian P (/x = 1) 
by the dashed line, and an optimized hermitian P 
(/j = 0.001) by the solid line. Broad shaded line is 
the exact analytic result. 

Closer inspection of Fig. 1, reveals that the sim- 
ulated expectation value of the Y quadrature, 
does not quite match the analytically predicted 
value for // = 1 for large times. This systematic 
error is due to non- vanishing boundary terms in 
the variable, making the change from master 
to Fokker-Planck equations inexact. The dis- 
crepancy can be reduced by using the optimized 
gauge with \x = 0.01, given by the solid line, al- 
though the sampling error increases. It is clear 
that further investigation into the trade-offs 
between reducing sampling error and reducing 
boundary term error is required. 



6 Final Comments 

The successful control and immense reduction 
of sampling error in the above one-mode exam- 
ple gives us confidence that the sampling error 
in the many-mode calculation can also be re- 



duced using this method, and BEC's can be sim- 
ulated after the point of condensation reached 
in [9]. The particular realization of the stochas- 
tic gauge idea discussed above is aimed toward 
the simulation of a BEC. However the approach 
is quite general, and may also be fruitful for 
simulations of many-mode higher dimensional 
bosonic systems. 
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